#include <stdio.h>
#include <string>
#include <sstream>
#include <iomanip>
#include "GridLoading\cgnsload.h"
#include "grid.h"
#include "gridgeneration.h"
#include "model.h"
#include "riemannsolvers.h"
#include "garbarukmodel.h"

template< typename T >
std::string int_to_hex( T i )
{
  std::stringstream stream;
  stream << "0x" 
         << std::setfill ('0') << std::setw(sizeof(T)*2) 
         << std::hex << i;
  return stream.str();
}

//Main program ))
int main(int argc, char *argv[]) {	

	//Load cgns grid
	Grid grid = LoadCGNSGrid("C:\\Users\\Erik\\Dropbox\\Science\\!Projects\\Grids\\FlatPlateMesh.cgns");	

	// Initialize medium model and place boundary conditions
	Model<Roe3DSolverPerfectGas, GarbarukTurbulentModel> model;

	//Set fluid properties
	model.SetGamma(1.4);
	model.SetCv(0.717e3);
	model.SetViscosity(1.7894e-05);

	//Set computational settings
	model.SetCFLNumber(0.3);

	////Bind computational grid
	model.BindGrid(grid);	

	////Initial conditions
	ConservativeVariables initValues(0);

	Vector velocity(0.1,0,0);
	double pressure = 1e5;
	double temperature = 300;

	initValues = model.PrimitiveToConservativeVariables(velocity, pressure, temperature, model.medium);
	//model.SetInitialConditions(initValues);
	//Determine plate start coordinate
	const double xPlateStart = 0.2;
	model.SetInitialConditionsBlasius(xPlateStart, velocity.x ,initValues.ro, initValues.roE);
	model.SaveToTechPlot("init.dat");

	//Boundary conditions
	//Inlet boundary
	Model<Roe3DSolverPerfectGas, GarbarukTurbulentModel>::InletBoundaryCondition InletBC(model);
	InletBC.setParams(1e5, 300, Vector(0.1,0,0));

	//Outlet boundary
	Model<Roe3DSolverPerfectGas, GarbarukTurbulentModel>::SubsonicOutletBoundaryCondition OutletBC(model);
	OutletBC.setParams(1e5);

	//Symmetry boundary
	Model<Roe3DSolverPerfectGas, GarbarukTurbulentModel>::SymmetryBoundaryCondition SymmetryBC(model);

	//No slip boundary
	Model<Roe3DSolverPerfectGas, GarbarukTurbulentModel>::NoSlipBoundaryCondition NoSlipBC(model);

	//Set boundary conditions
	/*model.SetBoundaryCondition("Left", InletBC);
	model.SetBoundaryCondition("Right", OutletBC);
	model.SetBoundaryCondition("Top", NoSlipBC);
	model.SetBoundaryCondition("Bottom", NoSlipBC);*/

	model.SetBoundaryCondition("INLET_3D", InletBC);
	model.SetBoundaryCondition("OUTLET_3D", OutletBC);
	model.SetBoundaryCondition("PLATE_3D", NoSlipBC);
	model.SetBoundaryCondition("SYMMETRY_3D", SymmetryBC);
	model.SetBoundaryCondition("TOP_LEFT_3D", OutletBC);
	model.SetBoundaryCondition("TOP_RIGHT_3D", OutletBC);

	//Set wall boundaries	
	model.SetWallBoundary("PLATE_3D", true);

	//Load solution
	std::string solutionFile = "solOutlet.txt";
	//model.LoadSolution("sol.txt");

	//Run simulation
	double totalTime = 0;	
	for (int i = 0; i < 100000000; i++) {
		model.Step();				
		std::cout<<"Interation = "<<i<<"\n";
		std::cout<<"TimeStep = "<<model.stepInfo.TimeStep<<"\n";
		std::cout<<"Residual = "<<model.stepInfo.Residual<<"\n";
		std::cout<<"TotalTime = "<<model.totalTime<<"\n";
		if (i % 1000 == 0) {
			model.SaveSolution(solutionFile);
			model.SaveToTechPlot("test.dat");
		};
		if (model.totalTime > 1e-2) break;
	};

	//Save result to techplot
	model.SaveSolution(solutionFile);
	model.SaveToTechPlot("test.dat");
}